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ABSTRACT 

'  When  describing  active- set  methods  for  linearly  constrained  optimisation,  it  is  often  con¬ 
venient  to  treat  all  constraints  in  a  uniform  manner.  However,  in  many  problems  the  linear  con¬ 
straints  include  simple  bounds  on  the  variables  as  well  as  general  constraints.  Special  treatment 
of  bound  constraints  in  the  implementation  of  an  active-set  method  yields  significant  advantages 
in  computational  effort  and  storage  requirements.  In  this  paper,  we  describe  how  to  perform  the 
constraint- related  steps  of  an  active-set  method  when  the  constraint  matrix  is  dense  and  bounds 
arc  treated  separately.  These  steps  involve  updates  to  the  TQ  factorisation  of  the  working  set  of 
constraints  and  the  Cholcsky  factorization  of  the  projected  Hessian  (or  Hessian  approximation). 
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1.  Introduction 

Constrained  optimization  problems  often  include  a  set  of  linear  inequality  constraints,  which  may 
be  written  in  several  different  forms.  We  consider  the  following  three: 

L Cl  Ax  >  6; 

LC2  Ax  —  b,  t  <  *  <  u; 

LC3 

For  convenience  we  shall  always  assume  that  A  is  a  matrix  with  m  rows  and  n  columns.  The 
dimensions  of  other  quantities  follow  in  each  case.  The  constraints  involving  A  arc  called  general 
constraints,  and  inequalities  of  the  form  t  <  x  <  u  are  called  simple  hounds  or  just  bounds. 
If  necessary,  some  of  the  components  of  l  or  u  may  be  taken  as  — oo  or  oo.  (Note  that  general 
equality  constraints  may  be  represented  in  LCl  by  extending  the  relations  to  include  equality, 
and  in  LC3  by  specifying  the  same  value  for  the  corresponding  elements  of  l  and  u.) 

The  forms  LCl  -  LC3  arc  equivalent,  in  the  sense  that  any  set  of  linear  inequality  constraints 
may  be  expressed  in  each  of  the  forms,  given  suitable  definition  of  A,  h,  l,  u  and  x.  The  primary 
feature  of  interest  in  LC2  is  that  general  inequality  constraints  are  converted  to  general  equality 
constraints  by  adding  slack  variables. 

The  most  popular  methods  for  treating  linear  inequality  constraints  are  called  active-set 
methods  (see  Section  2).  An  essential  characteristic  of  these  methods  is  that  they  maintain  a 
prediction  of  the  set  of  constraints  that  are  active  at  the  solution  (this  prediction  will  be  called 
the  working  set).  The  working  set  is  updated  by  adding  and  deleting  constraints  as  the  iterations 
proceed.  In  presenting  the  formal  description  of  an  active-set  method,  it  is  often  convenient 
to  treat  all  inequality  constraints  uniformly,  since  the  strategics  that  determine  changes  in  the 
working  set  are  usually  unaffected  by  whether  or  not  a  constraint  is  a  general  constraint  or  a 
simple  bound. 

Ill  any  implementation  of  an  active-set  method,  changes  in  the  working  set  involve  updates 
to  a  certain  matrix  C  associated  with  the  working  set  (and  often  to  other  matrices  as  well).  The 
data  structures  chosen  for  an  implementation  will  inevitably  bo  more  efficient  for  one  constraint 
form  than  another.  If  the  implementation  is  based  on  form  LCl,  changes  in  the  working  set  lead 
to  changes  in  the  rows  of  C,  while  with  form  LC2,  changes  in  the  working  set  lead  to  changes  in 
the  columns  of  C.  With  form  LC3,  the  changes  involve  both  the  rows  and  the  columns  of  C. 

Corresponding  changes  must  also  be  made  to  some  factorisation  of  C.  For  reasons  of 
simplicity,  past  implementors  have  dealt  almost  exclusively  with  forms  LCl  and  LC2.  Some 
examples  in  the  literature  follow. 

LCl,  dense  A:  Rosen  (1960)  and  many  others  for  nonlinear  programs  (see  Cill  and  Murray,  1974, 
for  further  references);  Stoer  (1971)  for  constrained  linear  least  squares. 

LCl,  sparse  A:  Buckley  (1975)  for  nonlinear  programs. 

I,C2,  dense  A:  Milllin  (1979)  for  constrained  linear  least  squares;  Bartels  (1980)  for  linear  pro¬ 
grams. 

LC2,  sparse  A:  Commercial  mathematical  programming  systems  for  linear  and  integer  programs; 
Murtagh  and  Saunders  (1978,  1982)  Tor  nonlinear  programs. 

The  disadvantage  of  implementations  based  on  LCl  or  LG2  arises  when  the  “natural”  state¬ 
ment  of  the  constraints  corresponds  most  closely  to  LC3  (i.e.,  when  there  is  a  significant  number 
of  bounds  in  L Cl  or  general  inequalities  in  LC2).  In  such  cases,  a  large  proportion  of  the  rows  or 
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columns  of  C  will  be  those  of  the  identity  matrix: 


respectively.  Maintaining  a  factorisation  of  the  whole  of  C  therefore  involves  more  than  the  ideal 
amount  of  storage  and  work.  (Certain  economics  do  arise  automatically  if  C  is  treated  as  a  sparse 
matrix,  but  much  of  the  objection  remains.) 

Implementations  based  on  LC3  effectively  take  advantage  of  the  above  structure  in  C.  Few 
authors  have  previously  considered  the  associated  complications.  Zoutendijk  (1070)  and  Powell 
(1975)  have  considered  how  changes  in  the  working  set  may  be  performed  when  C  is  square  (with 
varying  dimension)  and  its  inverse  is  updated  in  product  form.  Gill  and  Murray  (1973)  discussed 
the  nature  of  the  updates  required  in  a  non-sim plex  linear  programming  method  based  on  an 
orthogonal  factorization  of  C. 

In  this  paper  we  discuss  the  implementation  of  an  active-set  method  suited  to  constraints  of 
the  form  LC3,  with  A  treated  as  a  dense  matrix.  We  describe  how  to  update  the  TQ  factors  of 
the  matrix  C  and  the  Cliolcsky  factors  of  the  accompanying  projected  Hessian  (or  approximate 
Hessian).  The  procedures  have  been  implemented  in  computer  software  for  linear  and  quadratic 
programs  and  for  linearly  constrained  optimization,  as  described  in  Gill  et  a l.  (1982a, b).  The 
principal  advantages  in  dealing  with  LC3  are  as  follows. 

1.  The  matrix  to  be  factorized  has  dimension  mL  X  nFn,  where  mL  <  min{m,n)  and  nFH  <  n. 
Further,  mL  and  n,.„  arc  oflen  much  smaller  than  these  bounds. 

2.  When  finite  differences  are  used  to'  approximate  derivatives,  special  treatment  of  bounds  may 
lead  to  significant  economics  in  function  evaluations. 

3.  Certain  methods  for  semi-definite  and  indefinite  quadratic  programming  may  construct  a 
temporary  set  of  simple  bounds  in  order  to  begin  optimization.  For  such  methods,  the  ability 
to  handle  bounds  efficiently  is  crucial  even  if  the  original  problem  does  not  contain  bounds. 

The  first  advantage  is  best  illustrated  by  the  case  of  linear  programming.  Standard  implemen¬ 
tations  of  the  primal  simplex  method  (l)anlzig,  1963)  apply  to  constraints  in  the  form  LC. 2.  These 
arc  most  efficient  when  m  n,  since  the  matrix  to  be  factorized  is  always  m  X  m.  If  most  of  the  n 
variables  in  l,C2  arc  slack  variables,  the  standard  device  for  avoiding  gross  inefficiency  is  to  solve 
the  dual  problem.  In  contrast,  if  the  form  l,C3  is  assumed  when  implementing  the  simplex  method 
(the  most  famous  of  all  active-set  methods!),  then  maximum  efficiency  is  obtained  regardless  of 
the  ratio  of  m  to  n.  This  advantage  is  all  the  more  important  for  nonlinear  problems,  where  the 
device  of  solving  the  dual  is  not  necessarily  applicable  or  efficient. 

The  techniques  given  here  may  be  applied  to  active-set  methods  for  general  optimisation 
problems,  whenever  linear  and  nonlinear  constraints  are  treated  separately  —  particularly  in' 
methods  that  solve  a  sequence  of  quadratic  programming  subproblems  (o.g.,  Murray,  1969;  Diggs, 
1972;  Han,  1976;  Powell,  1977;  Murray  and  Wright,  1982)  or  linearly  constrained  subproblems 
(c.g.,  Rosen  and  Kreuser,  (972;  Robinson,  1972;  Murtagh  and  Saunders,  1982). 


2.  Overview  of  an  active- set  method 

Apart  from  the  requirement  of  feasibility,  the  optimality  conditions  for  a  constrained  problem 
involve  only  the  constraints  that  arc  active  (hold  with  equality)  at  the  solution.  Active-set  methods 
are  based  on  an  attempt  to  identify  the  constraints  that  arc  active  at  the  solution,  and  to  treat 
these  as  equality  constraints  in  the  subproblems  that  define  the  iterates.  The  temporary  equalities 
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are  used  to  reduce  the  dimensionality  of  the  minimisation.  In  a  typical  active-set  method,  the 
direction  of  search  is  computed  by  solving  a  (usually  simplified)  subproblem  in  which  a  subset  of 
the  problem  constraints  are  treated  as  5 qualities .  The  subset  of  the  problem  constraints  used  to 
compute  the  search  direction  will  be  called  the  working  set. 

Before  giving  a  detailed  description  of  the  special  treatment  of  bounds,  we  consider  some  of 
the  main  steps  in  an  active-set  method.  Our  concern  is  with  the  fc-th  iteration,  and  the  associated 
iterate  x*.  Wc  denote  by  Ck  the  matrix  whose  rows  are  the  constraints  in  the  current  working 
set,  by  tk  the  number  of  constraints  in  the  working  set  (the  number  of  rows  of  Ck),  and  by  gk  the 
gradient  of  the  function  to  be  minimised,  evaluated  at  z*.  The  matrix  Z*  is  defined  as  a  matrix 
whose  columns  span  the  null  space  of  C*  (i.e.,  C* Zk  =  0);  this  paper  is  primarily  concerned  with 
active-set  methods  in  which  Zk  is  stored  explicitly. 

The  major  operations  associated  with  the  current  working  set  are: 

(i)  formation  of  the  projected  gradient  Zjy^; 

(ii)  solution  of  the  linear  system 

%k ffk^kP*  =  —Zlt9k  (0 

for  the  (n  —  <n)-dimensional  vector  pM\ 

(iii)  calculation  of  the  search  direction  pk  =  Zkp x; 

(iv)  calculation  of  a  Lagrange  multiplier  estimate  X*  by  solving 

min  || Cl\k  -  VjkHl  (2) 

for  some  vector  vk. 

(These  quantities  may  he  computed  in  other  mathematically  equivalent  forms;  sec  Gill,  Murray 
and  Wright  (1981)  for  a  discussion  of  alternatives.) 

The  matrix  //*  in  (1)  usually  represents  second-derivative  information  about  the  objective 
function,  but  is  not  necessarily  stored  explicitly.  For  example,  fit  may  be  the  exact  Hessian  of 
the  objective  function  in  a  quadratic  program,  or  a  factorised  representation  of  the  Hessian  in  a 
linear  least-squares  problem.  In  some  methods,  Ilk  will  be  a  quasi-Newton  approximation  of  the 
Hessian  matrix,  or  Z'kIIkZ k  itself  will  be  approximated.  For  simplicity,  we  shall  always  refer  to 
Ilk  as  the  “Hessian”,  and  to  Z%IIkZk  as  the  “projected  Hessian”. 


3.  Representation  of  the  working  set  and  associated  factorizations 

Our  concern  in  this  section  is  with  the  factorisations  used  in  an  active-set  algorithm,  and  the 
elTcct  of  the  separate  treatment  of  bounds.  Wc  shall  assume  that  rank(C*)  =  t*,  i.e.  that  the 
rows  of  C*  arc  linearly  independent.  (In  practice,  this  condition  can  be  enforced  by  suitable 
choice  of  the  working  set;  see  Gill  et  a/.,  1982a.)  For  simplicity  of  notation,  wc  temporarily  drop 
the  subscript  k  associated  with  the  current  iteration. 

At  a  typical  iteration,  the  working  Xct  of  t  constraints  will  include  a  mixture  of  general 
constraints  and  hounds.  If  the  working  set  contains  any  simple  bounds,  those  variables  will  be 
fixed  on  the  corresponding  bounds  during  the  given  iteration;  all  other  variables  will  be  regarded 
as  free.  Wc  use  the  suffices  "FX”  and  “Fit”  to  denote  items  associated  with  the  two  types  of 
variable.  Suppose  that  C  contains  nrx  bounds  and  mL  general  constraints  (so  that  t  =  nrx  +  mt). 
Let  A  denote  the  matrix  whose  rows  are  the  mL  general  constraints  in  the  working  set,  and  let 
nrH  denote  the  number  of  free  variables  (nrR  =  n  —  nrx).  If  bounds  are  not  treated  separately, 
nrx  =  0,  npx  —  n,  and  mL  =  t. 
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In  the  implementation  of  an  active- set  method,  the  indices  of  the  free  variables  and  of 
the  general  constraints  in  the  working  set  may  be  stored  in  lists  (and  relevant  vectors  ordered 
accordingly).  Hence,  we  shall  assume  without  loss  of  generality  that  the  last  nrx  variables  are 
fixed.  The  matrix  of  constraints  in  the  working  set  can  then  be  written  as 

where  Af  „  is  an  mL  X  nFR  matrix,  and  IFX  denotes  an  nFX-dimensional  identity  matrix. 

The  first  matrix  that  must  be  available  in  order  to  perform  the  calculations  (i)  through  (iv) 
is  the  matrix  Z,  whose  columns  form  a  basis  for  the  set  of  vectors  orthogonal  to  the  rows  of  C. 
The  special  form  of  (3)  means  that  Z  also  has  a  special  form,  which  involves  only  the  columns  of 
A  corresponding  to  free  variables.  Let  nx  denote  n  —  t,  the  number  of  columns  of  Z.  An  n  X  n* 
matrix  Z  whose  columns  arc  orthogonal  to  the  rows  of  C  in  (3)  is  given  by 

*-(V>  ,4> 

where  ZFH  is  an  nrR  X  nx  matrix  whose  columns  form  a  basis  for  the  subspace  of  j4|rR  (i.e., 
^FR  =  0).  (If  mL  is  zero,  ZFR  is  the  nFn-dimcnsional  identity  matrix.) 

We  shall  obtain  the  matrix  Zm  in  (4)  from  a  variant  of  the  usual  orthogonal  factorisation 
which  we  shall  call  the  TQ  factorization.  (The  reasons  for  using  the  TQ  factorisation  will  be 
discussed  in  Section  5.3,  when  we  consider  procedures  for  updating  the  matrix  factorisations 
following  a  constraint  deletion.)  The  TQ  factorization  of  A,„  is  defined  by 


AmQ  —  (  0  T),  (5) 

where  Q  is  an  nFR  x  nFI,  orthonormal  matrix,  and  T  is  an  mL  X  rot  “reverse”  triangular  matrix 
such  that  Tij  =  0  for  *  +  j  <  mL.  We  illustrate  the  form  of  T  with  a  4  X  4  example: 

X 

T  = 

X  X 

V  x  x  x 

(Clearly,  T  is  simply  a  lower-triangular  matrix  with  its  columns  in  reverse  order.)  From  (5)  it 
follows  that  the  first  nz  columns  of  Q  can  be  taken  as  the  columns  of  the  matrix  ZrK.  We  denote 
the  remaining  columns  of  Q  by  Yrtt  (the  columns  of  YrR  form  an  orthogonal  basis  for  the  subepace 
of  vectors  spanned  by  the  rows  of  /lrH). 

The  TQ  factorization  for  the  full  matrix  C  (3)  has  the  form 

c(Q  °  W  0  n"  0  W°  0  'rx\  (6) 

\0  IrxJ  V  Am  Arx  /  V  0  0  /„  /  \  0  T  Atx  ) 

We  emphasize  that  the  usefulness  of  the  TQ  factorization  does  not  depend  on  separate 
treatment  of  bounds,  since  the  TQ  factorization  of  the  full  matrix  G  may  also  be  computed 
and  updated  using  the  procedures  to  be  described  in  an  implementation  based  on  the  form  LCl. 


5.  Implementation  and  storage 


6 


4.  Calculation  of  the  saarch  direction  and  Lagrange  multipliers 

The  calculations  in  Section  2  simplify  when  the  working  set  and  its  factorisation  have  the  special 
forms  (3),  (4)  and  (6).  From  (4)  it  follows  that  the  search  direction  has  the  form  p  =  (pJR  0)T. 
Further,  ZTg  =  ZjRgrR  and  ZTIIZ  =  Zj.^i rRZTR.  Consequently,  the  computation  of  pFR 
involves  three  steps:  forming  the  vector  Z^RgrR\  solving  the  linear  system 

^R^rK^rn  Pm  =  -ZrRgFR  (7) 

for  the  vector  pz;  and  forming  the  vector  prH  =  ZrRpM.  (In  certain  contexts,  such  aa  quadratic 
programming  and  linear  least-squares,  the  known  form  of  the  objective  function  allows  substantial 
savings  in  solving  (7).)  The  work  involved  is  reduced  as  nFX  increases;  therefore,  if  the  working  set 
contains  any  bound  constraints,  less  work  is  required  if  bound  constraints  arc  treated  separately. 

In  the  active-set  algorithms  of  interest,  the  matrix  in  (7)  is  assumed  to  be  positive  definite, 
and  thus  (7)  is  solved  using  the  Cliolesky  factorisation  of  ZFRII  t.RZFR  (sec,  c.g.,  Wilkinson,  1965; 
Stewart,  1973): 

^fR^FR^m  =  RTR,  (8) 

where  R  is  upper  triangular. 

Simplifications  also  arise  in  solving  equation  (2)  for  the  Lagrange  multiplier  estimates.  Let  X 
be  partitioned  into  an  m,. -vector  XL  (the  multiplier  estimates  corresponding  to  the  general  linear 
constraints)  and  an  nFX-vcctor  X,x  (the  multiplier  estimates  corresponding  to  the  active  bound 
constraints).  From  (6),  the  equations  defining  the  multipliers  are 


<?TX  = 


( 


0 

/rx 


aL  V  Vx \  =  (  aJr\l  \  /«Vh 
Alx  /  \  /  V  XFX  +  Ajx\L  J  V  vrx 


(9) 


where  as  means  “equal  in  the  least-squares  sense”. 

The  vector  \L  is  the  least-squares  solution  of  the  first  nF„  equations  of  (9)  (which  are 
compatible  if  ZFRvrn  —  0): 

ALK  **  »W 


It  follows  from  (6)  that  \L  may  be  obtained  by  forming  Yj„vt.R  and  solving  the  mL  X  mL  non- 
singular  reverse-triangular  system 


The  multiplier  estimates  .associated  with  the  bound  constraints  may  then  be  computed  directly 
by  substituting  in  the  remaining  equations  of  (9),  i.e. 


XFX  —  vFX  —  /FX\t. 


Note  that  if  mL  is  zero,  XFX  is  given  simply  by  vFX. 

The  number  of  multiplications  required  to  solve  (9)  in  the  manner  given  above  is  nmL  to 
form  YFRvrR  and  Ayx\lt  and  to  solve  the  reverse-triangular  system.  The  saving  in  work 
compared  to  treating  all  constraints  uniformly  is  £nFX  +  nrx(n  +  mL)  multiplications. 


B.  Implementation  and  storage 

When  implementing  an  active-set  method  based  on  the  TQ  and  Cholesky  factorizations,  the 
matrices  to  be  stored  include,  from  (6)  and  (8),  the  nFR  X  nFn  matrix  Q,  the  mL-dimensional 
reverse  triangular  matrix  T,  and  the  nF-dimensional  upper  triangular  matrix  It.  The  matrix  Q 
is  conceptually  partitioned  into  two  submatriccs  —  ZFR,  the  first  n*  columns,  and  F,-Hl  the  last 
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mL  columns,  i.e. 

Q  =  (Zrn  Y,n). 

Changes  in  the  working  set  will  cause  changes  in  these  four  matrices.  From  (6)  we  see  that  any 
transformations  applied  to  the  columns  of  yFn  will  also  be  applied  to  the  columns  of  T;  from 
(8)  it  follows  that  any  transformations  applied  to  the  columns  of  ZrR  will  also  be  applied  to  the 
columns  of  It.  (The  matrix  R  may  also  be  changed  in  other  ways  —  for  example,  by  a  low-rank 
modification  in  a  quasi-Newton  method.  However,  we  consider  only  the  effect  of  changes  in  the 
working  set.) 

In  the  implementation,  the  nrR  X  n,„  matrix  Q  is  stored  explicitly,  in  the  upper  left  corner 
of  sufficiently  large  array  ZY  (in  a  general  problem,  nr„  may  be  as  large  as  n).  The  dimensions 
of  It  and  T  are  complementary  (in  the  sense  that  n*  +  mL  —  n,„),  and  hence  both  matrices 
are  stored  in  the  upper  left  corner  of  a  single  array  RT.  The  matrix  R  is  stored  in  the  first  n, 
columns  (corresponding  to  the  columns  of  Zrn ),  and  the  matrix  T  in  the  following  mt  columns 
(corresponding  to  the  columns  of  With  this  storage  arrangement,  rotations  applied  to 

columns  of  ZY  can  be  applied  to  exactly  the  same  columns  of  RT.  We  have  chosen  to  store  the 
triangular  matrices  R  and  T  as  two-dimensional  arrays  (rather  than  in  a  “packed”  form),  so  that 
separate  subroutines  arc  not  required  to  apply  linear  algebraic  operations  to  triangular  matrices. 
(In  our  implementation  we  use  a  set  of  linear  algebra  subroutines  similar  to  the  BLAS  (Lawson 
et  ai.,  1979)  to  perform  these  operations.) 

The  following  diagram  illustrates  the  parallel  storage  arrangements  in  ZY  and  RT  for  the  case 
ni  =  aftd  mL  =  3.  The  elements  of  ZtR,  YFR,  R  and  T  are  denoted  by  z,  y,  r,  and  t, 
respectively. 


mL  nFX 


8.  Changes  in  the  working  set 

Unless  the  correct  active  set  is  known  a  priori,  the  working  set  must  be  modified  during  the 
execution  of  an  active-set  method,  by  adding  and  deleting  constraints.  Because  of  the  simple 
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nature  of  these  changes,  it  is  possible  to  update  the  necessary  matrix  factorizations  to  correspond 
with  the  new  working  set.  In  the  remainder  of  this  section,  we  consider  how  to  update  the  TQ 
factorization  (6)  and  the  Cholcsky  factorization  (8)  following  a  single  change  in  the  working  set. 
If  several  constraints  are  to  be  added  or  deleted,  the  procedures  arc  repeated  as  necessary. 

The  discussion  of  updates  will  assume  a  general  familiarity  with  the  properties  of  plane 
rotations.  Sequences  of  plane  rotations  are  used  to  introduce  zeros  into  appropriate  positions  of 
a  vector  or  matrix,  and  have  exceptional  properties  of  numerical  stability  (see,  c.g.,  Wilkinson, 
1965,  pp.  47-48). 

Wc  shall  illustrate  each  modification  process  using  sequences  of  simple  diagrams,  following 
the  conventions  of  Cox  (1981)  to  show  the  elTccts  of  the  plane  rotations.  Each  diagram  depicts 
the  changes  resulting  from  one  plane  rotation.  The  following  symbols  are  used: 
x  denotes  a  non-zero  clement  that  is  not  altered; 
m  denotes  a  non-zero  clement  that  is  modified; 
f  denotes  a  previously  zero  element  that  is  tilled  in; 

0  denotes  a  previously  non-zero  clement  that  is  annihilated;  and 
(or  blank)  denotes  a  zero  clement  that  is  unaltered. 

In  the  algebraic  representation  of  the  updates,  barred  quantities  will  represent  the  “new”  values. 


6.1.  Adding  a  general  constraint.  When  a  general  constraint  is  added  to  the  working  set,  its 
index  can  simply  be  placed  at  the  end  of  the  list  of  indices  or  general  constraints  in  the  working 
set.  Therefore,  without  loss  of  generality  wc  shall  assume  that  the  new  constraint  is  added  as 
the  last  row  of  A.  The  row  dimension  of  ArR  and  the  dimension  of  T  will  thus  increase  by  one, 
and  the  column  dimension  of  7,vn  will  decrease  by  one.  (Note  that  the  column  dimension  of  d,R 
is  unchanged.)  Let  a7  denote  the  new  row  of  A,  partitioned  into  (a^H  a7x).  Let  w1  denote  the 
vector  a^R Q,  and  partition  uiTaa  (to*  uix),  so  that  xti[  consists  of  the  first  n,  components  of  wT. 
From  (5),  it  follows  that 


ArnQ  ( £1  )Q  ~  _  «?)■ 


We  see  that  a  new  matrix  Q  can  be  obtained  by  applying  a  sequence  of  plane  rotations  on  the  right 
of  Q  to  transform  the  vector  iex  to  suitable  form;  the  transformed  matrix  Q  then  becomes  Q.  The 
sequence  of  rotations  Like  linear  combinations  of  the  elements  of  to  reduce  it  to  a  multiple 
(say,  y)  of  ej,  where  c„  denotes  the  n,-th  coordinate  vector.  The  rotations  arc  constructed  to 
alter  pairs  of  components  in  the  order  (1,2),  (2,3),  ...,  (rz,  —  l,ns),  as  indicated  in  the  following 
diagrams,  which  depict  the  vector  tej  as  it  is  reduced  to  ye^: 

( x  x  x  x  x )  -»  ( 0  m  x  x  x )  -*  (•  Omxx)  -*  (•  ■  Omx)  -»  (•  •  •  Om). 


The  effect  of  these  transformations  can  be  expressed  algebraically  as 


At  nQ 


T 

...T 


)  =  (»  n 


By  construction,  the  rotations  ill  P  affect  only  the  first  n,  columns  of  Q,  so  that  the  last  mt 
columns  of  Q  are  identical  to  those  of  Q,  and  the  first  n,  columns  of  Q  are  linear  combinations 
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of  the  first  n„  columns  of  Q.  Hence, 


—  (  %rn  y  )» 

where  y,  the  transformed  last  column  of  ZeHl  becomes  the  first  column  of  Y TR. 

The  plane  rotations  applied  to  ZFR  also  transform  the  Cholcsky  factor  R  of  the  projected 
Hessian.  The  chosen  order  of  the  rotations  in  P  means  that  each  successive  rotation  has  the 
elTcct  of  introducing  a  subdiagonal  clement  into  the  upper-triangular  matrix  R,  as  shown  in  the 
following  sequence  of  diagrams.  For  clarity,  we  again  show  the  vector  at  the  top  as  it  is 
reduced  to  (0  ~f)T\  the  matrices  underneath  represent  the  transformed  version  of  R. 
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Since  the  last  column  of  the  matrix  ZFRP  is  not  part  of  ZFR,  the  hist  column  of  RP  can  be 
discarded.  The  remaining  matrix  is  then  restored  to  upper-triangular  form  by  a  forward  sweep 
of  row  rotations  (say,  the  matrix  P),  which  is  applied  on  the  left  to  eliminate  the  subdiagonal 
elements,  as  shown  in  the  followir.g  diagrams. 

x  x  x  x  mmmm  x  x  x  x  x  x  x  x 

xxxx  Ommm  ■  m  m  m  -xxx 

xxx  — *  xxx  —*  0mm— ►  •  m  m 

xx  xx  xx  0m 

x  x  x  x 

Let  R  denote  the  matrix  RP  with  its  last  column  deleted;  then  we  have 

M-(o> 

where  R  is  upper-triangular.  Note  chat  the  rotations  in  P  affect  R,  but  not  Q  or  T. 

The  number  of  multiplications  associated  with  adding  a  general  constraint  includes  the 
following  (where  only  the  highest-order  term  is  given):  n*B  to  form  aFR Q;  3 n*  for  the  two  sweeps 
of  rotations  applied  to  R;  and  3nKBn*  to  transform  the  appropriate  columns  of  Q.  (We  assume 
the  three-multiplication  form  of  a  plane  rotation;  sec  Gill  et  a/.,  1 974.) 


6.2.  Adding  a  bound.  When  a  bound  constraint  is  added  to  the  working  set,  a  previously 
free  variable  becomes  lixed  on  its  bound.  Thus,  the  column  dimension  of  AFR,  the  column  and 
row  dimensions  of  ZF„  and  the  dimension  of  Q  arc  decreased  by  one.  The  dimension  of  T  is 
unaltered. 

We  assume  that  the  new  fixed  variable  corresponds  to  the  last  (nKB-th)  column  of  APR;  in 
practice,  the  index  of  the  variable  at  the  end  of  the  list  of  free  variables  is  moved  to  the  position 
of  the  newly  fixed  variable.  Let  eJR  denote  the  nrR-th  coordinate  vector.  Addition  of  the  np„-th 


i 
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variable  to  the  working  set  causes  the  vector  eJH  to  be  added  as  the  Brst  row  of  C,  i.e. 

=  (  (10) 

and,  in  effect,  moves  the  last  column  of  /tFR  into  AFX.  Let  wT  denote  the  nFR- th  row  of  Q,  and 
note  that  w  has  unit  Euclidean  length.  As  in  Section  6.1,  let  wT  be  partitioned  as  (iu^  u/^). 
After  adding  the  bound  to  the  working  set,  it  follows  from  (3)  and  (10)  that 


In  order  to  compute  the  updated  TQ  factorization,  the  first  row  of  the  matrix  on  the  right- 
hand  side  of  (1 1)  must  be  reduced  to  the  nFR-th  coordinate  vector.  This  is  achieved  by  a  forward 
sweep  of  plane  rotations  (say,  /’)  that  alter  columns  1  through  nFlt,  in  the  order  (1,2),  ..., 
(nFR  —  1,  nF„),  such  that  the  nFR-th  row  and  column  of  Q  are  transformed  to  the  nFB-th  coordinate 
vector.  The  effect  of  the  rotations  in  P  on  the  matrix  Q  can  be  represented  as 

QP  =  (  2fr  V™  )P  =  (  Q  °  )  =  (  ZrK  ?rn  °  Y 

V  o  i  /  V  o  o  i  J 

The  effect  of  the  rotations  in  P  on  R  and  T  is  best  understood  by  considering  them  in  two 

groups.  Firstly,  the  rotations  that  alter  columns  1  through  nz  of  Q  affect  the  columns  of  It 

exactly  as  described  in  Section  6.1,  and  a  set  of  row  rotations  are  then  applied  to  restore  the 
upper-triangular  form  of  R.  Secondly,  the  rotations  that  alter  columns  nz  through  nFR  of  Q 
cause  elements  to  be  added  above  the  reverse  diagonal  of  7',  as  shown  in  the  following  diagrams. 
The  vector  at  the  top  shows  the  order  of  the  rotations,  with  T  below. 
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The  first  mL  columns  of  the  transformed  T  become  T ,  and  the  remaining  column  becomes 
the  column  of  A  corresponding  to  the  new  fixed  variable. 

Let  I*x  denote  the  (nFX  +  l)-ditnensional  identity  matrix.  The  final  configuration  is  thus 
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as  desired. 

The  number  of  multiplications  associated  with  adding  a  bound  constraint  includes  all  those 
needed  to  add  a  general  constraint,  with  an  additional  to  modify  T. 


Optimisation  with  a  Mixture  of  Hound*  and  General  Constraint * 


6.3.  Deleting  a  general  constraint.  When  a  general  constraint  (say,  the  »-th)  is  deleted  from 
the  working  set,  the  row  dimension  of  Arn  and  the  dimension  of  T  arc  decreased  by  one,  and  the 
column  dimension  of  ZVR  is  increased  by  one.  (As  in  Section  6.1,  the  column  dimension  of  Arm 
is  not  altered.)  From  (5),  we  have 

ArnQ  —  (  0  S  ), 

where  5  is  an  (mt  —  1)  X  mL  matrix  such  that  rows  1  through  *  —  1  are  in  reverse  triangular 
form,  and  the  remaining  rows  have  one  extra  clement  above  the  reverse  diagonal.  In  order  to 
reduce  the  last  mf,  —  1  columns  of  S  to  the  desired  reverse  triangular  form  of  T,  a  backward 
sweep  of  plane  rotations  (say,  P)  is  applied  on  the  right,  to  take  linear  combinations  of  columns 
(mt  —  »'  +  1,  mL  —  *),  . . . ,  (2, 1).  The  effect  of  these  rotations  is  shown  in  the  following  diagrams 
for  the  case  mL  =  6  and  i  =  3: 
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This  transformation  may  be  expressed  as 

>Ue(  '  £  )=^r«0  =  (0  T). 


The  first  n*  columns  of  Q  are  not  affected  by  the  rotations  in  P,  and  hence  the  first  n,  columns 
of  Zrn  arc  identical  to  the  columns  of  ZFR.  The  matrix  ZFR  is  given  by 


=  {ZFR  Z  ),  (12) 

where  the  new  (last)  column  z  of  ZFH  is  a  linear  combination  of  the  relevant  columns  of  Yra- 
(When  *  =  mL,  no  reduction  at  all  is  needed,  and  z  is  just  the  first  column  of  YFR.) 

Because  of  the  form  (12),  the  new  projected  Hessian  matrix  ZlFnlI  TRZTK  is  given  by 

Zrnll.A.  =  tM »-(*?  ’)■  (l3> 

V  vT  0  j 

where  v  =  ZFRHrnz  and  0  =  zTIlFHz.  (Note  that  =  I{Fn.)  In  cases  when  HTR  or  RtR  is 
a  quasi-Newton  approximation,  the  new  row  and  column  of  (13)  may  need  to  be  approximated. 

Let  r  denote  the  solution  of  R?r  =  v;  if  2FRH rnZrn  is  positive  definite,  the  quantity  0  —  rTr 
must  be  positive.  In  this  case,  only  one  further  step  of  the  row-wise  Cholesky  factorisation  is 
needed  to  compute  the  new  Cholesky  factor  R,  which  is  given  by 


where  p2  =  0  —  rTr.  If  the  matrix  ZjK II TKZrR  is  not  positive  definite,  the  Cholesky  factorisation 
may  be  undefined  or  ill-conditioned,  and  other  techniques  should  be  used  to  modify  the  factorisa¬ 
tion  without  excessive  additional  computation  or  loss  of  numerical  stability  (e.g.,  see  Gill  et  a I., 
1082a,  for  techniques  applicable  to  quadratic  programming). 
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The  number  of  multiplications  associated  with  deleting  the  t-th  general  constraint  includes 
the  following  (where  only  the  highest-ordcr  term  is  given):  *)*  to  operate  on  T j  3nr„(mi,— i) 

to  transform  Q;  n\R  to  form  IIrRz;  nFRn,  to  form  ZjHHFlxz;  and  to  compute  the  additional 
row  of  the  Cholcsky  factor.  It  is  clearly  advantageous  to  delete  constraints  at  the  end  of  the  list 
of  general  constraints  in  the  working  set;  hence,  the  indices  of  general  equality  constraints  are 
always  placed  at  the  beginning  of  the  list. 

The  justification  for  using  the  TQ  factorisation  arises  from  this  part  of  an  active-set  method. 
From  a  theoretical  viewpoint,  ZFR  would  remain  an  orthogonal  basis  for  the  null  space  of  AFR 
regardless  of  the  position  in  which  the  new  column  appeared.  However,  in  order  to  update  the 
Cholesky  factors  efficiently,  the  new  column  must  appear  after  the  columns  of  ZrH  (otherwise, 
(13)  would  not  hold).  The  TQ  factorisation  has  an  implementation  advantage  because  the  new 
column  of  ZFH  automatically  appears  in  the  correct  position  after  deletion  or  a  constraint  from 
the  working  set.  With  other  alternatives,  the  housekeeping  associated  with  the  update  of  ft  is 
more  complicated.  For  example,  in  an  implementation  based  on  the  LQ  factorisation,  the  new 
column  might  be  moved  to  the  end  of  ZFR,  or  a  list  could  be  maintained  of  the  locations  of  the 
columns  of  ZFR;  another  alternative  is  to  store  the  columns  of  ZFR  in  reverse  order  (see  Gill  and 
Murray,  1977). 


6.4.  Deleting  a  bound.  When  a  bound  constraint  is  deleted  from  the  working  set,  a  previously 
fixed  variable  becomes  free.  In  this  case,  the  column  dimension  of  •/!,„,  the  column  and  row 
dimensions  of  ZFR  and  the  dimension  of  Q  -Arc  increased  by  one;  the  dimension  of  T  remains 
unaltered.  In  practice,  the  index  of  the  newly  freed  variable  is  added  at  the  end  of  the  list  of 
free  variables.  Hence,  we  assume  without  loss  of  generality  that  the  (nrR  +  l)-th  variable  is  to  be 
freed  from  its  bound,  so  that  0  is  defined  by  deleting  row  nFR  +  1  of  C.  In  effect,  the  boundary 
between  ylKR  and  Arx  is  “shifted”  by  one  column;  this,  corresponds  to  augmenting  Q  by  a  row 
and  column  of  the  identity,  and  reducing  by  one  the  dimension  of  the  identity  matrix  associated 
with  the  fixed  variables. 

Let  a  denote  the  column  of  A  corresponding  to  the  newly  freed  variable,  and  let  l~x  denote 
the  (nFX  -  l)-dimcnsional  identity  matrix.  The  result  of  deleting  the  bound  constraint  is  then 

o(Q  0  )  =  (  0  0  '")f»  °  °  W0  0  0 

\0  IFX  J  \  ArR  a  Arx  J  \  Q  0  J-^)  \0  T  a  AFX 

To  reduce  the  augmented  matrix  ( T  a)  to  the  desired  form  (0  T),  a  backward  sweep  of 

column  plane  rotations  (say,  P)  is  applied  in  the  order  (mt  +  l,mt),. .. ,  (2,  l),  as  shown  in  the. 
following  diagrams: 
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The  rotations  in  P  affect  columns  n,  +  1  through  .nFR  +  1  of  the  augmented  Q.  The  first  n, 
columns  of  ZrH  are  thus  simply  those  of  ZFR,  with  a  row  of  icros  added  at  the  bottom.  It  follows 
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that  ZTR  is  given  by 


(14) 


where  the  last  element  of  z  will  be  nonxero. 

The  effect  of  freeing  the  (nKR  +  l)-th  variable  is  to  augment  HFR  by  the  row  and  column 
corresponding  to  the  newly  released  variable.  Since  (14)  is  similar  to  (12),  the  Cholcsky  factors  of 
the  new  projected  Hessian  can  be  obtained  from  the  existing  factors  by  performing  one  further  step 
of  the  factorization  as  before  (assuming  that  the  updated  projected  Hessian  is  positive  definite). 

The  number  of  multiplications  associated  with  deleting  a  bound  constraint  includes  the 
following,  where  only  the  highest-order  term  is  given:  to  operate  on  T;  3 mLnrR  to  transform 

Q>  n?«i  to  form  f/FRz;  nFRnx  to  form  Z^RHrRz;  and  to  update  R. 


References 

Bartels,  R.  H.  (1980).  A  penalty  linear  programming  method  using  reduced-gradient  basis- 
exchange  techniques,  Linear  Algebra  and  its  Applies.  20,  pp.  17-32. 

Biggs,  M.  C.  (1972).  “Constrained  minimization  using  recursive  equality  quadratic  programming”, 
in  Numerical  Methods  for  Non-Linear  Optimisation  (F.  A.  Lootsma,  ed.),  pp.  411-428, 
Academic  Press,  London  and  New  York. 

Buckley,  A.  G.  (1975).  An  alternative  implementation  of  Goldfarb’s  minimization  algorithm, 
Math.  Prog.  8,  pp.  207-231. 

Cox,  M.  G.  (1981).  The  least  squares  solution  of  overdetermined  linear  equations  having  band  or 
augmented  band  structure,  IMA  J.  Numcr.  Anal.  1,  pp.  3  -22. 

Dantzig,  G.  B.  (1963).  Linear  programming  and  extensions,  Princeton  University  Press,  Princeton, 
New  Jersey. 

Gill,  P.  E.,  Golub,  G.  II.,  Murray,  W.,  and  Saunders,  M.  A.  (1974).  Methods  for  updating  matrix 
factorizations,  Math.  Comp.  28,  pp.  505-535. 

Gill,  P.  E.  and  Murray,  W.  (1973).  A  numerically  stable  form  of  the  simplex  algorithm,  Lin.  Alg. 
Applies.  7,  pp.  99  138. 

Gill,  P.  E.  and  Murray,  W.  (1974).  Newton-type  methods  for  unconstrained  and  linearly  con¬ 
strained  optimization,  Math.  Prog.  28,  pp.  311-350. 

Gill,  P.  E.,  and  Murray,  W.  (1977).  “Linearly  constrained  problems,  including  linear  and  quad¬ 
ratic  programming,”  in  The  State  of  the  Art  in  Numerical  Analysis  (D.  Jacobs,  ed.),  pp. 
313-363,  Academic  Press,  London  and  New  York. 

Gill,  P.  E.,  Murray,  W.,  Saunders,  M.  A.,  and  Wright,  M.  H.  (1982a).  The  design  and  implemen¬ 
tation  of  a  quadratic  programming  Algorithm,  Report  SOL  82-8,  Department  of  Operations 
Research,  Stanford  University,  California. 

Gill,  P.  E.,  Murray,  W.,  Saunders,  M.  A.,  and  Wright,  M.  II.  (1982b).  Software  for  linearly 
constrained  optimization,  Report  SOI,  82-9,  Department  of  Operations  Research,  Stanford 
University,  California. 

Gill,  P.  E.,  Murray,  W.  and  Wright,  M.  II.  (1981).  Practical  Optimisation,  Academic  Press, 
London  and  New  York. 


References 


13 


Han,  S.-P.  (1976).  Superiinearly  convergent  variable  metric  algorithms  for  general  nonlinear 
programming  problems,  Math.  Prog.  11,  pp.  263-282. 

Lawson,  C.  L.,  Hanson,  R.  J.,  Kincaid,  D.  R.,  and  Krogh,  F.  T.  (1979).  Basic  linear  algebra 
subprograms  for  Fortran  usage,  ACM  Trans.  Math.  Software  6,  pp.  308-325. 

Mifflin,  R.  (1979).  A  stable  method  for  solving  certain  constrained  least-squares  problems,  Math. 
Prog.  16,  pp.  141-158. 

Murray,  W.  (1969).  “An  algorithm  for  constrained  minimization”,  in  Optimization  (R.  Fletcher, 
ed.),  pp.  247-258,  Academic  Press,  London  and  New  York. 

Murray,  W.  and  Wright,  M.  H.  (1982).  Computation  of  the  search  direction  in  constrained 
optimisation  algorithms,  Math.  Prog.  Study  16,  pp.  62-83. 

Murtagh,  B.  A.  and  Saunders,  M.  A.  (1978).  Large-scale  linearly  constrained  optimization,  Math. 
Prog.  14,  pp.  41-72. 

Murtagh,  B.  A.  and  Saunders,  M.  A.  (1982).  A  projected  Lagrangian  algorithm  and  its  implemen¬ 
tation  for  sparse  nonlinear  constraints,  Math.  Prog.  Study  1C,  pp.  84-117. 

Powell,  M.  J.  D.  (1977).  A  fast  algorithm  for  nonlinearly  constrained  optimization  calculations, 
Report  DAMTP  77/NA  2,  University  of  Cambridge,  England. 

Powell,  S.  (1975).  A  development  of  the  product  form  algorithm  for  the  simplex  method  using 
reduced  transformation  vectors,  Math.  Prog.  Study  4,  pp.  93-107. 

Robinson,  S.  M.  (1972).  A  quadratically  convergent  algorithm  for  general  nonlinear  programming 
problems,  Math.  Prog.  3,  pp.  145-156. 

Rosen,  J.  B.  (1960).  The  gradient  projection  method  for  nonlinear  programming,  Part  I  —  linear 
constraints,  SIAM  J.  Appl.  Math.  8,  pp.  181-217. 

Rosen,  J.  B.  and  Kreuser,  J.  (1972).  “A  gradient  projection  algorithm  for  nonlinear  constraints”, 
in  Numerical  Methods  for  Non-Linear  Optimization  (F.  A.  Loolsma,  ed.),  pp.  297  300, 
Academic  Press,  London  and  New  York. 

Stewart,  G.  W.  (1973).  Introduction  to  Matrix  Computations,  Academic  Press,  London  and  New 
York. 

Stocr,  J.  (1971).  On  the  numerical  solution  of  constrained  least-squares  problems,  SIAM  J.  N timer. 
Anal.  8,  pp.  382-411. 

Wilkinson,  J.  H.  (1965).  The  Algebraic  Eigenvalue  Problem,  Oxford  University  Press. 

Zoutendijk,  G.  (1970).  “A  product-form  algorithm  using  contracted  transformation  vectors”,  in 
Integer  and  Nonlinear  Programming  (J.  Abadie,  ed.),  pp.  511-523,  North-IIolland,  Amster¬ 
dam. 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  FAOK  (Wkam  SW.SWllS. 


1  REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

A  RECIPIENT’S  CATALOG  NUMBER 

«.  TITLE  (andkabUUa)  1 

Procedures  for  Optimization  Problems  with  a 
Mixture  of  Bounds  and  General  Linear  Con¬ 
straints. 

Technical  Report 

S.  PERFORMING  ORO.  REPORT  NUMBER 

t.  author?*; 

Philip  E.  Gill,  Walter  Murray 

Michael  A.  Saunders,  Margaret  H.  Wright 

i.  dONTRACT  OR  GRANT  NUMBER?*; 

N00014-75-C-0267 

DAAG29-81-K-0156 

S.  PERFORMING  ORGANIZATION  NAME  ANO  ADDRESS 

Department  of  Operations  Research  -  SOL 

Stanford  University 

Stanford,  CA  94305 

II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Office  of  Naval  Research  -  Dept,  of  the  Navy 

800  N.  Quincy  Street 

Arlington,  VA  22217 

It.  REPORT  date 

May,  1982 

IS.  NUMBER  OF  PAGES 

13  PP, 

U.  MONITORING  AGENCY  NAME  *  ADDRESS?!/  dtttatent  ham  CanUatHnt  Olllea) 

U.S.  Army  Research  Office 

P.0.  Box  12211 

Research  Triangle  Park,  NC  27709 

IS.  SECURITY  CLASS,  (at  thla  report) 

UNCLASSIFIED 

1  16.  DISTRIBUTION  STATEMENT  (ol  tfif •  Rmport) 

This  document  has  been  approved  for  public  release  and  sale; 
Its  distribution  Is  unlimited. 


IT.  DISTRIBUTION  STATEMENT  (a!  the  abetted  entered  In  Black  39,  II  dlllarant  Bam  Raman) 


•s.  supplementary  notes 


IS.  KEY  WORDS  (Canlhtua  an  rarataa  aide  II  necaaaaay  and  Identity  by  black  tnmrbar) 

-Inear  constraints;  active-set  methods;  updating  matrix  factorizations; 
computer  software;  orthogonal  factorization;  TQ  factorization;  Cholesky 
factorization. 


20.  ABSTRACT  (ConfimM  on  reverN  tide  If  n«e««Mry  «d  Identify  by  block  wbir) 

See  other  side. 


jjCMWTY  CLAMWCATIOH  OF  THU  WAOIfWhi  Ml  Mnt*n4) 


SOL  82-6:  PROCEDURES  FOR  OPTIMIZATION  PROBLEMS  WITH  A  MIXTURE  OF 
BOUNDS  AND  GENERAL  LINEAR  CONSTRAINTS  (May  1982,  13  pp) 

When  describing  active-set  methods  for  linearly  constrained 
optimization,  it  is  often  convenient  to  treat  all  constraints  in  a 
uniform  manner.  However,  in  many  problems  the  linear  constraints 
include  simple  bounds  on  the  variables  as  well  as  general 
constraints.  Special  treatment  of  bound  constraints  in  the 
implementation  of  an  active- set  method  yields  significant  advantages 
in  computational  effort  and  storage  requirements.  In  this  paper,  we 
describe  how  to  perform  the  constraint-related  steps  of  an  active-set 
method  when  the  constraint  matrix  is  dense  and  bounds  are  treated 
separately.  These  steps  involve  updates  to  the  TQ  factorization  of 
the  working  set  of  constraints  and  the  Cholesky  factorization  of  the 
projected  Hessian  (or  Hessian  approximation). 


